/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/* edp.c -- Direct evaluation of membrane density profiles.
   Profiles are described by histograms along the lines of rdf histograms.
   Memory is allocated in module startUp.c, function AllocArrays().
   References: RDF algorithms in Allen-Tildesley, 6.2 and Rapaport, 4.3 */

#include "dataStructs.h"

extern const VecR *siteCoords; 
extern const Site *site; 
extern const VecR region; 
extern const int nSites, stepCount, sizeHistEdp, limitEdp;
extern const real soluteMass, deltaT, timeNow, edpDeltaZ; // bin size along the bilayer normal (z-axis) computed in step.c
extern int countEdp, halfSizeHistEdp; // globals, defined & initialised in startUp.c
extern real meanEdpSlabVol, meanEdpDeltaZ, meanEdpBoxHeight; // globals, defined and initialised in startUp.c

/* global variables */
real *soluteHistEdp, *watHistEdp, *cholHistEdp, *phosHistEdp, *glyHistEdp, *estHistEdp, *totalHistEdp;
real *CH2_3HistEdp, *CHCH2_2HistEdp, *CH3CH2_2HistEdp;  

void ZeroEdpHistograms(); // external function declaration

static void PrintEdp()
{
  real distanceFromBilayerCentre;
  int bin;

  FILE *edpFilePtr;
  char edpSnapFile[ 20 ];
  char time_ns[ 9 ]; // string containing timeNow, e.g. "100", in ns

  sprintf( time_ns, "%d", ( int ) ( 0.5 + timeNow *  TIME_ns ) );
  edpSnapFile[ 0 ] = '\0';
  strcat( edpSnapFile, "edp-" );
  strcat( edpSnapFile, time_ns );
  strcat( edpSnapFile, ".dat" );
  edpFilePtr = fopen( edpSnapFile, "w" );
  
  for ( bin = 0; bin < sizeHistEdp; bin++ ) {
    distanceFromBilayerCentre = ( bin + 0.5 ) * meanEdpDeltaZ - meanEdpBoxHeight / 2.; 
    fprintf( edpFilePtr, "%7.3f", distanceFromBilayerCentre ); /* [nm] */
    fprintf( edpFilePtr, "%12.3f", totalHistEdp[ bin ] ); /* [electrons/nm^3] */
    fprintf( edpFilePtr, "%12.3f", watHistEdp[ bin ]  );
    fprintf( edpFilePtr, "%12.3f", cholHistEdp[ bin ] );
    fprintf( edpFilePtr, "%12.3f", phosHistEdp[ bin ] );
    fprintf( edpFilePtr, "%12.3f", glyHistEdp[ bin ] );
    fprintf( edpFilePtr, "%12.3f", estHistEdp[ bin ] );
    fprintf( edpFilePtr, "%12.3f", CH2_3HistEdp[ bin ] );
    fprintf( edpFilePtr, "%12.3f", CHCH2_2HistEdp[ bin ] );
    fprintf( edpFilePtr, "%12.3f", CH3CH2_2HistEdp[ bin ] );
    /*     fprintf( edpFilePtr, "%12.3f", soluteHistEdp[ bin ] ); */
    fprintf( edpFilePtr, "\n" );
  }
  fclose( edpFilePtr );
}

void EvalEdp()
{
  real slabVolume,
    zDistCom; // site COM z-distance
  int n, bin; // "bin" is the histogram bin number

  /* assigning electrons - total is 374 for DMPC (Tristram-Nagle, 2002), 438 for DSPC, 434 for DOPC (Liu & Nagle, PRE, 2004) */
  int neH2O = 10; /* number of electrons per water molecule */
  int neH = 1; /* number of electrons of H atom */
  int neC = 6; /* number of electrons of C atom */
  int neN = 7;
  int neO = 8;
  int neP = 15; 
  int neCH = 7;
  int neCH2 = 8; 
  int neCH3 = 9;
  int neAmine = 3*neH + neN + 2*neCH2; /* total number of electrons of amine site (DOPE) */
  int neCholine = 3*neCH3 + 2*neCH2 + neN; /* total number of electrons of choline site (DOPC) */
  int nePhosphate = neP + 4*neO; /* total number of electrons of phosphate site */
  int neGlycerol = 2*neCH2 + neCH + 2*neO; /* total number of electrons of glycerol site */
  int neEster = 2*neCH2 + neC + neO; /* total number of electrons of ester site */
  int neCH2_3    = 3*neCH2; 
  int neCHCH2_2  = 2*neCH2 + neCH; 
  int neCH3CH2_2 = 2*neCH2 + neCH3; 

  slabVolume = region.x * region.y * edpDeltaZ;   /* computing the volume over which the density is worked out */
  
  DO_SITE { // computing edp histograms: sorting sites in bins from 0 to sizeHistEdp - 1 
    zDistCom = siteCoords[ n ].z; 
    bin = zDistCom / edpDeltaZ + halfSizeHistEdp;    // printf ("***site[ %d ].type = %d: bin = %d\n", n, site[ n ].type, bin);
    switch ( site[ n ].lipidUnit ) {
    case WATER:  
      watHistEdp[ bin ] += neH2O;  
      break;
    case CHOLINE_OR_AMINE: 
      if ( site[n].type == CHOLINE_TYPE ) 
	cholHistEdp[ bin ] += neCholine; /* DOPC choline headgroup */
      else if ( site[n].type == AMINE_TYPE ) 
	cholHistEdp[ bin ] += neAmine; /* DOPE amine headgroup */
      break;
    case PHOSPHATE: 
      phosHistEdp[ bin ] += nePhosphate; 
      break;
    case GLYCEROL: 
      glyHistEdp[ bin ] += neGlycerol; 
      break;
    case ESTER_A: 
    case ESTER_B: 
      estHistEdp[ bin ] += neEster; 
      break;
    case TAIL_A1:
    case TAIL_A4:
    case TAIL_B1:
    case TAIL_B4:
      CH2_3HistEdp[ bin ] += neCH2_3;
      break;
    case TAIL_A2:
    case TAIL_A3:
    case TAIL_B2:
    case TAIL_B3:
      CHCH2_2HistEdp[ bin ] += neCHCH2_2;
      break;
    case TAIL_A5:
    case TAIL_B5:
      CH3CH2_2HistEdp[ bin ] += neCH3CH2_2;
      break;
    default: printf("WARNING: default reached in EvalEdp switching.\n"); 
    }
  }

  meanEdpSlabVol += slabVolume; // accumulating - strictly necessary only at constant pressure to account for volume changes
  meanEdpDeltaZ += edpDeltaZ;  // edpDeltaZ computed in SingleStep()
  meanEdpBoxHeight += region.z; 
    
  ++countEdp;

  if ( limitEdp == countEdp ) {  // if happy with the amount of data collected
    for ( bin = 0; bin < sizeHistEdp; bin++ ) { // normalising to have an averaged value
      /*       soluteHistEdp[ bin ]   /= meanEdpSlabVol;  */
      watHistEdp[ bin ]       /= meanEdpSlabVol; 
      cholHistEdp[ bin ]      /= meanEdpSlabVol; 
      phosHistEdp[ bin ]      /= meanEdpSlabVol; 
      glyHistEdp[ bin ]      /= meanEdpSlabVol; 
      estHistEdp[ bin ]       /= meanEdpSlabVol; 
      CH2_3HistEdp[ bin ]    /= meanEdpSlabVol; 
      CHCH2_2HistEdp[ bin ]  /= meanEdpSlabVol; 
      CH3CH2_2HistEdp[ bin ] /= meanEdpSlabVol; 
      totalHistEdp[ bin ] = watHistEdp[ bin ] + cholHistEdp[ bin ] + phosHistEdp[ bin ] /* accumulating total */
	+ glyHistEdp[ bin ] + estHistEdp[ bin ] + CH2_3HistEdp[ bin ] + CHCH2_2HistEdp[ bin ] + CH3CH2_2HistEdp[ bin ]; 
    }
    
    meanEdpDeltaZ /= countEdp;
    meanEdpBoxHeight /= countEdp; 
    
    PrintEdp();
    countEdp = 0;
    meanEdpSlabVol = 0.;
    ZeroEdpHistograms();
  }
}

